*
* $Id$
*
* $Log: prepre.F,v $
* Revision 1.1.1.1  2002/06/16 15:18:43  hristov
* Separate distribution  of Geant3
*
* Revision 1.1.1.1  1999/05/18 15:55:22  fca
* AliRoot sources
*
* Revision 1.1.1.1  1995/10/24 10:22:03  cernlib
* Geant
*
*
#include "geant321/pilot.h"
*CMZ :  3.21/02 28/03/94  01.30.59  by  S.Giani
*-- Author :
#if defined(CERNLIB_HPUX)
$OPTIMIZE OFF
#endif
*$ CREATE PREPRE.FOR
*COPY PREPRE
*
*=== prepre ===========================================================*
*
      SUBROUTINE PREPRE ( WEE, NNEUT, NPROT, NHOLE, ICYCL )
 
#include "geant321/dblprc.inc"
#include "geant321/dimpar.inc"
#include "geant321/iounit.inc"
*
*----------------------------------------------------------------------*
*----------------------------------------------------------------------*
*
#include "geant321/balanc.inc"
#include "geant321/eva0.inc"
#include "geant321/fheavy.inc"
#include "geant321/finuc.inc"
#include "geant321/nucdat.inc"
#include "geant321/nucgeo.inc"
#include "geant321/parevt.inc"
#include "geant321/parnuc.inc"
#include "geant321/part.inc"
#include "geant321/resnuc.inc"
*
      REAL RNDM(2)
      COMMON / FKCOSP / C1ST (3), C2ND (3), LEMISS
      LOGICAL LEMISS
      COMMON / FKCMCY / NPCYCL (20,2), IEVT, LOUT
      COMMON / FKPLOC / IABCOU
*
      DIMENSION ERNCM  (2), EPSMX  (2), AMMAFT (2), AMNAFT (2),
     &          ZAFTR  (2), EEXMNM (2), EEXDLT (2), EEXANW (2),
     &          EPSEEX (2), EPSDLT (2), EPSANY (2), EPSFIX (29),
     &          C (3), UMULEG (0:2), ACOLEG (0:2), NPART (2),
     &          NPRFIX (2)
      LOGICAL LSPREJ, LNUCTS
*
      NPART (1) = NNEUT
      NPART (2) = NPROT
      IF ( LHLFIX ) THEN
         NHLFIX = NHLEXP
         NHOLE  = NHOLE - NHLFIX
         EHLFIX = 0.D+00
         DO 50 IHOLE = 1, NHLFIX
            EHLFIX = EHLFIX + HOLEXP (IHOLE)
   50    CONTINUE
         IF ( NHOLE .GT. 0 ) THEN
            RHOIMP = ( RHOIMP * ( NHOLE + NHLFIX ) - RHOEXP ) / NHOLE
            EKFIMP = ( EKFIMP * ( NHOLE + NHLFIX ) - EKFEXP ) / NHOLE
         ELSE
            RHOIMP = RHOAVE
            EKFIMP = EKFAVE
         END IF
      ELSE
         NHLFIX = 0
         EHLFIX = 0.D+00
      END IF
      NHINI  = NHOLE
      IF ( PTRES .LT. ANGLGB ) THEN
         UMO2  = ERES * ERES
         UMO   = SQRT (UMO2)
         EKRES = 0.D+00
         GAMCM = 1.D+00
         ETAX  = 0.D+00
         ETAY  = 0.D+00
         ETAZ  = 0.D+00
         PXORI = 0.D+00
         PYORI = 0.D+00
         PZORI = 0.D+00
         PCMORI = 0.D+00
         CALL RACO ( CXAXCM, CYAXCM, CZAXCM )
      ELSE
         UMO2  = ( ERES - PTRES ) * ( ERES + PTRES )
         UMO   = SQRT (UMO2)
         EKRES = ERES - UMO
         GAMCM = ERES  / UMO
         ETACM = PTRES / UMO
         ETAX  = PXRES / UMO
         ETAY  = PYRES / UMO
         ETAZ  = PZRES / UMO
         ETAPCM = ETAX * PXORI + ETAY * PYORI + ETAZ * PZORI
         PHELP = EKORI + AM (KPORI) - ETAPCM / ( GAMCM + 1.D+00 )
         PCMSX = PXORI - ETAX * PHELP
         PCMSY = PYORI - ETAY * PHELP
         PCMSZ = PZORI - ETAZ * PHELP
         PCMORI = SQRT ( PCMSX**2 + PCMSY**2 + PCMSZ**2 )
         CXAXCM = PCMSX / PCMORI
         CYAXCM = PCMSY / PCMORI
         CZAXCM = PCMSZ / PCMORI
      END IF
      BBTAR = IBTAR
      ZZTAR = ICHTAR
      TVCMS  = UMO - AMNRES
      EHLFIX = MIN ( TVCMS, EHLFIX )
      TVEFF  = TVCMS - EHLFIX
      IF ( TVEFF .LE. 0.D+00 ) GO TO 7000
      TEMNU = SQRT ( TVCMS / ANOW * ALEVEL  )
      NEXMX = SQRT ( 2.D+00 * ANOW * TVCMS / ALEVEL )
      NPTOT  = NPART (1) + NPART (2)
      IF ( NHOLE + NPTOT + NHLFIX .GE. NEXMX .OR.
     &     NPTOT .GT. NINT (0.5D+00 * ANOW) ) GO TO 7000
      ANPTOT = NPTOT
      ANPROT = NPART (2)
      ANNEUT = NPART (1)
      AVEBIN = ( ( BBTAR - ZZTAR - ACOLL + ZCOLL ) * BNENRG (2)
     &       + ( ZZTAR - ZCOLL ) * BNENRG (1) ) / ( BBTAR - ACOLL )
      IAAFT = IBRES-1
      IZAFT = ICRES
      AAFTR     = IAAFT
      ZAFTR (1) = IZAFT
      AMMAFT(1) = AAFTR * AMUAMU + 1.D-03 * FKENER (AAFTR,ZAFTR(1))
      AMNAFT(1) = AMMAFT (1) - ZAFTR (1) * AMELEC + ELBNDE (IZAFT)
      CALL EEXLVL ( IAAFT, IZAFT, EEXDLT (1), EEXMNM (1), EEXDUM )
      BNDGAV = BNENRG (2)
      EPSMX  (1) = UMO - AMNRES - BNDGAV - EHLFIX
      EPSDLT (1) = UMO - AMNRES - BNDGAV - EEXDLT (1)
      EPSEEX (1) = UMO - AMNRES - BNDGAV - EEXMNM (1)
      IZAFT = ICRES - 1
      ZAFTR (2) = IZAFT
      AMMAFT(2) = AAFTR * AMUAMU + 1.D-03 * FKENER (AAFTR,ZAFTR(2))
      AMNAFT(2) = AMMAFT (2) - ZAFTR (2) * AMELEC + ELBNDE (IZAFT)
      CALL EEXLVL ( IAAFT, IZAFT, EEXDLT (2), EEXMNM (2), EEXDUM )
      BNDGAV = BNENRG (1)
      EPSMX  (2) = UMO - AMNRES - BNDGAV - EHLFIX
      IF ( EPSMX (1) + EPSMX (2) .LT. 2.D+00 * TEMNU ) GO TO 7000
      EPSDLT (2) = UMO - AMNRES - BNDGAV - EEXDLT (2)
      EPSEEX (2) = UMO - AMNRES - BNDGAV - EEXMNM (2)
      IF ( NP .LE. NP0 ) THEN
         IF ( KPORI .EQ. 8 ) THEN
            EEXANW (1) = EEXDLT (1)
            EPSANY (1) = MIN ( EPSDLT (1), EPSMX (1) )
            EEXANW (2) = 0.D+00
            EPSANY (2) = EPSMX  (2)
         ELSE IF ( KPORI .EQ. 1 ) THEN
            EEXANW (1) = 0.D+00
            EPSANY (1) = EPSMX  (1)
            EEXANW (2) = EEXDLT (2)
            EPSANY (2) = MIN ( EPSDLT (2), EPSMX (2) )
         ELSE
            EEXANW (1) = 0.D+00
            EPSANY (1) = EPSMX (1)
            EEXANW (2) = 0.D+00
            EPSANY (2) = EPSMX (2)
         END IF
      ELSE
         EEXANW (1) = 0.D+00
         EPSANY (1) = EPSMX (1)
         EEXANW (2) = 0.D+00
         EPSANY (2) = EPSMX (2)
      END IF
      SIGIN0 = PI * ( R0SIGK * RMASS (IBRES-1) )**2
      NEMISS = 0
      MNUCTS = 0
 1000 CONTINUE
         NPRFIX (1) = 0
         NPRFIX (2) = 0
         JEMFIX = 0
         EPRFIX = 0.D+00
         IF ( .NOT. LEMISS ) ICYCL = ICYCL + 1
         IF ( NPTOT .LE. 0 ) GO TO 4600
         IF ( NNUCTS .GT. MNUCTS ) THEN
            MNUCTS = MNUCTS + 1
            JNUCTS = INUCTS (MNUCTS)
            JEMIS0 = 2 - ABS (KPNUCL(JNUCTS)) / 8
            JEMIS1 = JEMIS0
            JDEMIS = 1
            EKFSAV = EKFIMP
            RHOSAV = RHOIMP
            TMPRHO = 0.001D+00 * RHOAVE
            RHOIMP = MAX ( RHNUCL (JNUCTS), TMPRHO )
            TMPEKF = 0.001D+00 * EKFAVE
            EKFIMP = MAX ( EKFNUC (JNUCTS), TMPEKF )
            LNUCTS = .TRUE.
         ELSE
            LNUCTS = .FALSE.
            CALL GRNDM(RNDM,1)
            IF ( RNDM (1) .LT. ANNEUT / ANPTOT ) THEN
               JEMIS0 = 1
               JEMIS1 = 2
               JDEMIS = 1
            ELSE
               JEMIS0 = 2
               JEMIS1 = 1
               JDEMIS = -1
            END IF
         END IF
         NEXTOT = NPTOT + NHOLE
         RHORED = ACOLL / BBTAR
         IF ( LGDHPR ) THEN
            EKFCNN = ( ZNOW * EKFCEN (1) + ( ANOW - ZNOW ) * EKFCEN (2)
     &             ) / ANOW
            RHOWEI = 0.5D+00
            IF ( LNUCTS ) THEN
               AEFREF = EKFIMP
               AEFRAV = EKFIMP
               RHONUC = RHOIMP
            ELSE IF ( NPTOT .LE. 2 .AND. ICYCL .EQ. 1 ) THEN
               AHLTOT = MAX ( NHOLE, 1 )
               WEIGH1 = NHINI / AHLTOT
               WEIGH2 = 1.D+00 - WEIGH1
               AEFREF = WEIGH1 * EKFIMP + WEIGH2 * EKFAVE
               RHONUC = WEIGH1 * RHOIMP + WEIGH2 * RHOAVE
               AEFRAV = WEIGH1 * EKFIMP + WEIGH2 * EKFAVE
               RHONUC = RHORED * ( ( AHLTOT + 1.D+00 - RHOWEI )
     &                * RHONUC + RHOWEI * RHOAVE ) / ( AHLTOT + 1.D+00 )
               AEFRAV = ( AHLTOT * AEFRAV + EKFAVE )
     &                / ( AHLTOT + 1.D+00 )
            ELSE
               AHLTOT = MAX ( NHOLE, 1 )
               WEIGH1 = NHINI / AHLTOT
               WEIGH2 = 1.D+00 - WEIGH1
               AEFREF = WEIGH1 * EKFIMP + WEIGH2 * EKFAVE
               RHONUC = WEIGH1 * RHOIMP + WEIGH2 * RHOAVE
               AEFRAV = WEIGH1 * EKFIMP + WEIGH2 * EKFAVE
               RHONUC = RHORED * ( ( AHLTOT + 1.D+00 - RHOWEI )
     &                * RHONUC + RHOWEI * RHOAVE ) / ( AHLTOT + 1.D+00 )
               AEFRAV = ( AHLTOT * AEFRAV + EKFAVE )
     &                / ( AHLTOT + 1.D+00 )
            END IF
            CLAMDI = 1.0D+00
         ELSE
            EKFCNN = AEFRMX
            AEFREF = 0.5D+00 * AEFRMX
            AEFRAV = AEFREF
            RHONUC = 0.5D+00**1.5D+00 * RHONU0
            CLAMDI = 0.5D+00
         END IF
         DO 4500 JEMISS = JEMIS0, JEMIS1, JDEMIS
            IF ( NPART (JEMISS) .LE. 0 ) GO TO 4500
            IEMISS = 3 - JEMISS
            RATEC0 = 2.D+00 * AMNUCL (IEMISS) / PLABRC**3 / PISQ
            BNDGEN = AMNAFT (JEMISS) + AMNUCL (IEMISS) - AMNRES
            IF ( .NOT. LNUCTS .AND. EPRFIX .GT. ANGLGB ) THEN
               TVEFF  = TVEFF - EPRFIX
               NPTOT  = NPTOT - NPRFIX (1) - NPRFIX (2)
               NPART (1) = NPART (1) - NPRFIX (1)
               NPART (2) = NPART (2) - NPRFIX (2)
               JEMFIX = JEMISS
               ERNCM0 = ERNCM (JEMISS)
               EPSMX0 = EPSMX (JEMISS)
               BNDGAV = BNENRG (IEMISS)
               EPSMX (JEMISS) = UMO - AMNRES - BNDGAV - EHLFIX - EPRFIX
            ELSE
               JEMFIX = 0
            END IF
            BNDGAV = BNENRG (IEMISS)
            EPSMIN = BNDGEN - BNDGAV
            BNDHLP = BNENRG (IEMISS)
            IF ( EPSMX (JEMISS) - EPSMIN .LT. TEMNU ) GO TO 4500
            IF ( EHLFIX + EPRFIX .LE. 0.5D+00 * EEXANW (JEMISS) ) THEN
               UUMIN  = TVEFF - 0.5D+00 * ( EPSANY (JEMISS)
     &                + EPSMX (JEMISS) ) - BNDGAV
            ELSE
               UUMIN  = TVEFF - EPSMX (JEMISS) - BNDGAV
            END IF
            UUMAX  = TVEFF - BNDGEN
            BNDUSE = BNENRG (IEMISS)
            IF ( LNUCTS ) THEN
               EKNNUC = ENNUC  (JNUCTS) - AMNUCL (IEMISS) + BNDGEN
     &                - RSTNUC (JNUCTS)
               EPSHLP = BNDHLP - BNDGEN + EPSMIN
               EKNNUC = EKNNUC + EPSHLP
               IF ( EKNNUC .LE. EPSMIN ) THEN
                  GO TO 4500
               ELSE IF ( EKNNUC .LE. 0.D+00 ) THEN
                  DELTAE = 0.5D+00 * ( - EPSMIN - EKNNUC )
                  EKNNUC = EKNNUC + DELTAE
               ELSE
                  DELTAE = 0.D+00
               END IF
               ENNUC (JNUCTS) = EKNNUC + AMNUCL (IEMISS)
               PNUCCO = SQRT ( EKNNUC * ( EKNNUC + 2.D+00
     &                * AMNUCL (IEMISS) ) )
               PHELP = PNUCCO / PNUCL (JNUCTS)
               PCMSX = PXNUCL (JNUCTS) * PHELP
               PCMSY = PYNUCL (JNUCTS) * PHELP
               PCMSZ = PZNUCL (JNUCTS) * PHELP
               ETAPCM = ETAX * PCMSX + ETAY * PCMSY + ETAZ * PCMSZ
               PHELP = ENNUC (JNUCTS) - ETAPCM / ( GAMCM + 1.D+00 )
               PCMSX = PCMSX - ETAX * PHELP
               PCMSY = PCMSY - ETAY * PHELP
               PCMSZ = PCMSZ - ETAZ * PHELP
               PCMST = SQRT ( PCMSX**2 + PCMSY**2 + PCMSZ**2 )
               EPS = GAMCM * ENNUC (JNUCTS) - ETAPCM - AMNUCL (IEMISS)
               IF ( EPS .LE. EPSMIN ) THEN
                  GO TO 4500
               ELSE IF ( EPS .LT. ANGLGB * ENNUC (JNUCTS) ) THEN
                  GO TO 4500
               END IF
               EPS   = EPS - DELTAE
               C (1) = PCMSX / PCMST
               C (2) = PCMSY / PCMST
               C (3) = PCMSZ / PCMST
               EPS = MIN ( EPS, HLFHLF * ( EPSMX (JEMISS)
     &             + EPSANY (JEMISS) ) )
               IF ( JEMISS .EQ. 2 ) THEN
                  FLKCOU = DOST ( 1, ZAFTR (JEMISS) )
                  CCOUL  = DOST ( 3, ZAFTR (JEMISS) )
                  ETHRES = FLKCOU * COULBH * ( ZNOW - 1.D+00 )
     &                   / RMASS (IBRES-1)
               END IF
               FLAMMX = 1.D+00
               NTENT  = 1
               ISAMPL = 0
               LSPREJ = .FALSE.
               GO TO 3000
            END IF
            IF ( UUMIN .GE. UUMAX ) GO TO 4500
            EPSINV = 0.5D+00 * ( UMO2 - ( AMNUCL (IEMISS)
     &             + AMNAFT (JEMISS) )**2 ) / AMNAFT (JEMISS)
            EPSWLL = EPSINV + AEFRAV + BNDUSE + AMNUCL (IEMISS)
            BETWLL = SQRT ( 1.D+00 - AMNUSQ (IEMISS) / EPSWLL**2 )
            EKEWLL = EPSWLL - AMNUCL (IEMISS)
            IF ( JEMISS .EQ. 1 ) THEN
               AALPHA = 0.76D+00 + 2.2D+00 / RMASS (IBRES-1)
               BBETAA = ( 2.22D+00 / RMASS (IBRES-1)**2 - 0.05D+00 )
     &                / AALPHA * 1.D-03
               TMP07  = 0.07D+00
               EPSEFF = MIN ( EPSINV, TMP07 )
               SIGINV = 0.1D+00 * XINNEU ( EPSEFF, ZAFTR (JEMISS),
     &                  BBETAA )
               SFTYMX = 1.1D+00
               GSNGLV = 1.5D+00 * ( ANOW - ZNOW ) / EKFCEN (IEMISS)
               IF ( EKEWLL .LE. 0.04D+00 ) THEN
                  SIGMAP = 3.D+03 * PI / ( 1.206D+03 * EKEWLL + (
     &                   -1.86D+00 + 0.09415D+03 * EKEWLL
     &                   + 0.0001306D+06 * EKEWLL**2 )**2 ) + 1.D+03
     &                   * PI / ( 1.206D+03 * EKEWLL + ( 0.4223D+00
     &                   + 0.13D+03 * EKEWLL )**2 )
                  IF ( EKEWLL .LT. 0.02D+00 ) THEN
                     SIGMAN = 0.3333333333333333D+00 * SIGMAP
                  ELSE
                     SIGMAN = 10.63D+00 / BETWLL**2 - 29.92D+00 / BETWLL
     &                      + 42.9D+00
                  END IF
               ELSE
                  SIGMAP = 34.10D+00 / BETWLL**2 - 82.2D+00 / BETWLL
     &                   + 82.2D+00
                  SIGMAN = 10.63D+00 / BETWLL**2 - 29.92D+00 / BETWLL
     &                   + 42.9D+00
               END IF
               SIGMAP = 0.1D+00 * SIGMAP
               SIGMAN = 0.1D+00 * SIGMAN
               SGNUNU = ( ( ANOW - ZNOW - ANNEUT ) * SIGMAN
     &                + ( ZNOW - ANPROT ) * SIGMAP )
     &                / ( ANOW - ANPTOT )
            ELSE
               EPSCOU = UMO - AMNAFT (JEMISS) - AMNUCL (IEMISS)
               FLKCOU = DOST ( 1, ZAFTR (JEMISS) )
               CCOUL  = DOST ( 3, ZAFTR (JEMISS) )
               ETHRES = FLKCOU * COULBH * ( ZNOW - 1.D+00 )
     &                / RMASS (IBRES-1)
               TMP07  = 0.07D+00
               EPSEFF = MIN ( EPSINV, TMP07 )
               SIGINV = 1.D-01 * XINPRO ( EPSEFF, ZAFTR (JEMISS),
     &                                    ETHRES )
               SFTYMX = 1.2D+00
               GSNGLV = 1.5D+00 * ZNOW / EKFCEN (IEMISS)
               IF ( EKEWLL .LE. 0.04D+00 ) THEN
                  SIGMAN = 3.D+03 * PI / ( 1.206D+03 * EKEWLL + (
     &                   -1.86D+00 + 0.09415D+03 * EKEWLL
     &                   + 0.0001306D+06 * EKEWLL**2 )**2 ) + 1.D+03
     &                   * PI / ( 1.206D+03 * EKEWLL + ( 0.4223D+00
     &                   + 0.13D+03 * EKEWLL )**2 )
                  IF ( EKEWLL .LT. 0.02D+00 ) THEN
                     SIGMAP = 0.3333333333333333D+00 * SIGMAN
                  ELSE
                     SIGMAP = 10.63D+00 / BETWLL**2 - 29.92D+00 / BETWLL
     &                      + 42.9D+00
                  END IF
               ELSE
                  SIGMAN = 34.10D+00 / BETWLL**2 - 82.2D+00 / BETWLL
     &                   + 82.2D+00
                  SIGMAP = 10.63D+00 / BETWLL**2 - 29.92D+00 / BETWLL
     &                   + 42.9D+00
               END IF
               SIGMAP = 0.1D+00 * SIGMAP
               SIGMAN = 0.1D+00 * SIGMAN
               SGNUNU = ( ( ANOW - ZNOW - ANNEUT ) * SIGMAN
     &                + ( ZNOW - ANPROT ) * SIGMAP )
     &                / ( ANOW - ANPTOT )
            END IF
            ZITA   = AEFRAV / EKEWLL
            IF ( ZITA .LE. 0.5D+00 ) THEN
               PZITA = 1.D+00 - 1.4D+00 * ZITA
            ELSE
               PZITA = 1.D+00 - 1.4D+00 * ZITA + 0.4D+00 * ZITA
     &               * ( 2.D+00 - 1.D+00 / ZITA )**2.5D+00
            END IF
            RATEC  = SIGINV * RATEC0 * EPSINV / GSNGLV
            IF ( RATEC .LE. 0.D+00 ) GO TO 4500
            ALAMDC = BETWLL / RATEC
            ALAMDI = 1.D+00 / ( SGNUNU * RHONUC * PZITA ) + MAX (TWOTWO
     &             * PI * PLABRC * BETWLL * EPSWLL / AMNUSQ (IEMISS),
     &             RZNUCL )
            AMUTOT = 1.D+00 / ALAMDI + 1.D+00 / ALAMDC
            RATEI  = BETWLL / ALAMDI
            FLAMMX = SFTYMX * RATEC / ( RATEC + CLAMDI * RATEI )
            IF ( FLAMMX .LT. 0.D+00 ) THEN
               GO TO 4500
            END IF
            LSPREJ = .NOT. ( NHOLE .GT. 2 .OR. NPTOT .GT. 4 )
            IF ( LSPREJ ) THEN
               IF ( NPTOT .EQ. 1 ) THEN
                  IF ( NHOLE .EQ. 0 ) THEN
                     EPS = EPSMX (JEMISS)
                     FLAMMX = 1.D+00
                     NTENT  = 1
                     ISAMPL = 0
                     LSPREJ = .FALSE.
                     GO TO 3000
                  ELSE IF ( NHOLE .EQ. 1 ) THEN
                     ISAMPL = 1
                  ELSE IF ( NHOLE .EQ. 2 ) THEN
                     ISAMPL = 2
                  ELSE
                     ISAMPL = 7
                  END IF
               ELSE IF ( NPTOT .EQ. 2 ) THEN
                  IF ( NHOLE .EQ. 0 ) THEN
                     ISAMPL = 7
                  ELSE IF ( NHOLE .EQ. 1 ) THEN
                     ISAMPL = 3
                  ELSE IF ( NHOLE .EQ. 2 ) THEN
                     ISAMPL = 4
                  ELSE
                     ISAMPL = 7
                  END IF
               ELSE IF ( NPTOT .EQ. 3 ) THEN
                  IF ( NHOLE .EQ. 0 ) THEN
                     ISAMPL = 7
                  ELSE IF ( NHOLE .EQ. 1 ) THEN
                     ISAMPL = 5
                  ELSE IF ( NHOLE .EQ. 2 ) THEN
                     ISAMPL = 6
                  ELSE
                     ISAMPL = 7
                  END IF
               ELSE
                  ISAMPL = 7
               END IF
            ELSE
               ISAMPL = 7
            END IF
            GO TO ( 2100, 2200, 2300, 2400, 2500, 2600, 2700 ), ISAMPL
 2100       CONTINUE
               LSPREJ = .FALSE.
               AEFRPA = AEFREF
               UUMAX  = MIN ( AEFRPA, UUMAX )
               AITOT  = ( AEFRPA**3 - ( AEFRPA + UUMIN - UUMAX )**3 )
     &                / AEFRPA**3
               IF ( UUMAX .GT. TVEFF ) THEN
                  AIPRO = MIN ( AITOT, ONEONE ) * FLAMMX
                  AITOT = AITOT * FLAMMX
               ELSE
                  AITOT = AITOT * FLAMMX
                  AIPRO = AITOT
               END IF
            GO TO 2800
 2200       CONTINUE
               LSPREJ = .FALSE.
               AEFRPA = AEFREF
               UUMAX  = MIN ( TWOTWO * AEFRPA, UUMAX )
               IF ( TVEFF .LE. AEFRPA ) THEN
                  EDENOM = TVEFF**2 * 0.5D+00
               ELSE IF ( TVEFF .LE. 2.D+00 * AEFRPA ) THEN
                  EDENOM = ( TVEFF**2 - 2.D+00 * ( TVEFF - AEFRPA )**2 )
     &                   * 0.5D+00
               ELSE
                  EDENOM = AEFRPA**2
               END IF
               UUMN2 = UUMIN**2
               IF ( UUMAX .GT. AEFRPA ) THEN
                  UUDIV = AEFRPA
                  UUDV2 = UUDIV**2
                  AIPR1 = 0.5D+00 * UUDV2
                  AIPR2 = 0.5D+00 * ( UUMAX - UUDIV ) * ( 3.D+00 * UUDIV
     &                  - UUMAX )
                  IF ( UUMIN .LE. UUDIV ) THEN
                     AITO1 = 0.5D+00 * ( UUDV2 - UUMN2 )
                     AITO2 = AIPR2
                     AIHLP = 0.D+00
                  ELSE
                     AITO1 = 0.D+00
                     AIHLP = 0.5D+00 * ( UUMIN - UUDIV ) * ( 3.D+00
     &                     * UUDIV - UUMIN )
                     AITO2 = AIPR2 - AIHLP
                  END IF
               ELSE
                  UUDIV = UUMAX
                  UUDV2 = UUDIV**2
                  AIPR1 = 0.5D+00 * UUDV2
                  AITO1 = 0.5D+00 * ( UUDV2 - UUMN2 )
                  AITO2 = 0.D+00
                  AIPR2 = 0.D+00
               END IF
               AITOT = AITO1 + AITO2
               AITO1 = AITO1 / AITOT
               AITO2 = AITO2 / AITOT
               AITOT = AITOT * FLAMMX / EDENOM
               IF ( UUMAX .GT. TVEFF ) THEN
                  AIPRO = MIN ( ( AIPR1 + AIPR2 ) / EDENOM, ONEONE )
     &                  * FLAMMX
               ELSE
                  AIPRO = ( AIPR1 + AIPR2 ) * FLAMMX / EDENOM
               END IF
            GO TO 2800
 2300       CONTINUE
               LSPREJ = .FALSE.
               AEFRPA = AEFREF
               IF ( TVEFF .LE. AEFRPA ) THEN
                  EDENOM = 0.25D+00 * TVEFF * TVEFF
               ELSE
                  EDENOM = 0.25D+00 * AEFRPA * ( 2.D+00 * TVEFF
     &                   - AEFRPA )
               END IF
               FCHLP = 0.25D+00 * NPART (JEMISS)
               IF ( UUMAX .LE. AEFRPA ) THEN
                  UUDIV = UUMAX
                  AIPR1 = FCHLP * UUDIV**2
                  AITO1 = AIPR1 - FCHLP * UUMIN**2
                  AITO2 = 0.D+00
                  AIPR2 = 0.D+00
               ELSE
                  UUDIV = AEFRPA
                  AIPR1 = FCHLP * UUDIV**2
                  AIPR2 = FCHLP * ( UUMAX - UUDIV ) * UUDIV
                  IF ( UUMIN .GT. UUDIV ) THEN
                     AITO1 = 0.D+00
                     AIHLP = FCHLP * ( UUMIN - UUDIV ) * UUDIV
                     AITO2 = AIPR2 - AIHLP
                  ELSE
                     AIHLP = 0.D+00
                     AITO1 = AIPR1 - FCHLP * UUMIN**2
                     AITO2 = AIPR2
                  END IF
               END IF
               AITOT = AITO1 + AITO2
               AITO1 = AITO1 / AITOT
               AITO2 = AITO2 / AITOT
               IF ( UUMAX .GT. TVEFF ) THEN
                  DDNPAR = NPART(JEMISS)
                  AIPRO = MIN ( ( AIPR1 + AIPR2 ) / EDENOM,
     &                    DDNPAR ) * FLAMMX
               ELSE
                  AIPRO = ( AIPR1 + AIPR2 ) * FLAMMX / EDENOM
               END IF
            GO TO 2800
 2400       CONTINUE
               LSPREJ = .TRUE.
               AEFRPA = AEFREF
               IF ( TVEFF .LE. AEFRPA ) THEN
                  LSPREJ = .FALSE.
                  FEMAX  = 1.D+00
                  ISAMPL = 7
               ELSE IF ( TVEFF .LE. 2.D+00 * AEFRPA ) THEN
                  FEMAX = TVEFF**3 / ( TVEFF**3 - 2.D+00 * ( TVEFF
     &                  - AEFRPA )**3 )
               ELSE
                  FEMAX = 0.1666666666666667D+00 * TVEFF**3
     &                  / ( AEFRPA**2 * ( TVEFF - AEFRPA ) )
               END IF
               UUMNR  = ( UUMIN / TVEFF )**(NEXTOT-1)
               UUMXR  = ( UUMAX / TVEFF )**(NEXTOT-1)
               AITOT  = NPART (JEMISS) * FLAMMX * ( UUMXR - UUMNR )
     &                * FEMAX
               AIPRO  = NPART (JEMISS) * FLAMMX * UUMXR * FEMAX
            GO TO 2800
 2500       CONTINUE
               LSPREJ = .FALSE.
               AEFRPA = AEFREF
               IF ( TVEFF .LE. AEFRPA ) THEN
                  EDENOM = TVEFF**3 / 36.D+00
               ELSE
                  EDENOM = ( TVEFF**3 - ( TVEFF - AEFRPA )**3 )
     &                   / 26.D+00
               END IF
               IF ( UUMAX .GT. AEFRPA ) THEN
                  UUDIV = AEFRPA
                  UUDV3 = UUDIV**3
                  UUMN3 = UUMIN**3
                  FCHLP = NPART(JEMISS) / 36.D+00
                  AIPR1 = FCHLP * UUDV3
                  AIPR2 = 3.D+00 * FCHLP * UUMAX * UUDIV * ( UUMAX
     &                  - UUDIV )
                  IF ( UUMIN .LE. UUDIV ) THEN
                     AITO1 = AIPR1 - FCHLP * UUMN3
                     AITO2 = AIPR2
                     AIHLP = 0.D+00
                  ELSE
                     AITO1 = 0.D+00
                     AIHLP = 3.D+00 * FCHLP * UUMIN * UUDIV * ( UUMIN
     &                     - UUDIV )
                     AITO2 = AIPR2 - AIHLP
                  END IF
               ELSE
                  UUDIV = UUMAX
                  UUDV3 = UUDIV**3
                  UUMN3 = UUMIN**3
                  FCHLP = NPART(JEMISS) / 36.D+00
                  AIPR1 = FCHLP * UUDV3
                  AITO1 = AIPR1 - FCHLP * UUMN3
                  AIPR2 = 0.D+00
                  AITO2 = 0.D+00
               END IF
               AITOT = AITO1 + AITO2
               AITO1 = AITO1 / AITOT
               AITO2 = AITO2 / AITOT
               AITOT = AITOT * FLAMMX / EDENOM
               IF ( UUMAX .GT. TVEFF ) THEN
                  DDNPAR = NPART(JEMISS)
                  AIPRO = MIN ( ( AIPR1 + AIPR2 ) / EDENOM,
     &                    DDNPAR ) * FLAMMX
               ELSE
                  AIPRO = ( AIPR1 + AIPR2 ) * FLAMMX / EDENOM
               END IF
            GO TO 2800
 2600       CONTINUE
               AEFRPA = AEFREF
               IF ( TVEFF .LE. AEFRPA ) THEN
                  EDENOM = TVEFF**4 / 288.D+00
               ELSE IF ( TVEFF .LE. 2.D+00 * AEFRPA ) THEN
                  EDENOM = ( TVEFF**4 - 2.D+00 * ( TVEFF - AEFRPA )**4 )
     &                   / 288.D+00
               ELSE
                  EDENOM = AEFRPA**2 * ( 0.5D+00 * TVEFF**2 - TVEFF
     &                   * AEFRPA + 7.D+00 / 12.D+00 * AEFRPA**2 )
     &                   / 12.D+00
               END IF
               FCHLP = NPART (JEMISS) / 288.D+00
               IF ( UUMAX .GT. 2.D+00 * AEFRPA ) THEN
                  LSPREJ = .TRUE.
                  UUDIV = AEFRPA
                  UUDI2 = 2.D+00 * AEFRPA
                  UUDV4 = UUDIV**4
                  UUMN4 = UUMIN**4
                  UUD24 = UUDI2**4
                  AIPR1 = FCHLP * UUDV4
                  AIPR2 = FCHLP * ( UUD24 - UUDV4 )
                  AIPR3 = 12.D+00 * FCHLP * AEFRPA**2 * UUMAX
     &                  * ( UUMAX - UUDI2 )
                  IF ( UUMIN .LE. UUDIV ) THEN
                     AITO1 = AIPR1 - FCHLP * UUMN4
                     AITO2 = AIPR2
                     AITO3 = AIPR3
                     AIHLP = 0.D+00
                     AIHL2 = 0.D+00
                  ELSE IF ( UUMIN .LE. UUDI2 ) THEN
                     AITO1 = 0.D+00
                     AIHLP = FCHLP * ( UUMN4 - UUDV4 )
                     AITO2 = AIPR2 - AIHLP
                     AITO3 = AIPR3
                     AIHL2 = 0.D+00
                  ELSE
                     AITO1 = 0.D+00
                     AITO2 = 0.D+00
                     AIHL2 = 12.D+00 * FCHLP * AEFRPA**2 * UUMIN
     &                     * ( UUMIN - UUDI2 )
                     AITO3 = AIPR3 - AIHL2
                  END IF
               ELSE IF ( UUMAX .GT. AEFRPA ) THEN
                  LSPREJ = .TRUE.
                  UUDIV = AEFRPA
                  UUDI2 = UUMAX
                  UUDV4 = UUDIV**4
                  UUMN4 = UUMIN**4
                  UUD24 = UUDI2**4
                  AIPR1 = FCHLP * UUDV4
                  AIPR2 = FCHLP * ( UUD24 - UUDV4 )
                  IF ( UUMIN .LE. UUDIV ) THEN
                     AITO1 = AIPR1 - FCHLP * UUMN4
                     AITO2 = AIPR2
                     AIHLP = 0.D+00
                  ELSE
                     AITO1 = 0.D+00
                     AIHLP = FCHLP * ( UUMN4 - UUDV4 )
                     AITO2 = AIPR2 - AIHLP
                  END IF
                  AIPR3 = 0.D+00
                  AITO3 = 0.D+00
               ELSE
                  LSPREJ = .FALSE.
                  UUDI2 = UUMAX
                  UUDIV = UUMAX
                  UUDV4 = UUDIV**4
                  UUMN4 = UUMIN**4
                  AIPR1 = FCHLP * UUDV4
                  AITO1 = AIPR1 - FCHLP * UUMN4
                  AITO2 = 0.D+00
                  AIPR2 = 0.D+00
                  AITO3 = 0.D+00
                  AIPR3 = 0.D+00
               END IF
               AITOT = AITO1 + AITO2 + AITO3
               AITO1 = AITO1 / AITOT
               AITO2 = AITO2 / AITOT
               AITO3 = AITO3 / AITOT
               AITOT = AITOT * FLAMMX / EDENOM
               AIPRO = ( AIPR1 + AIPR2 + AIPR3 ) * FLAMMX / EDENOM
            GO TO 2800
 2700       CONTINUE
               LSPREJ = .FALSE.
               TVLEFF = LOG(TVEFF)
               IF(UUMAX.LE.0.) THEN
                  UULMXR=-100
               ELSE
                  UULMXR = (LOG(UUMAX)-TVLEFF)*(NEXTOT-1)
               ENDIF
               IF(UULMXR.LT.-88) THEN
                  UUMNR = 0.
                  UUMXR = 0.
               ELSE
                  UUMXR = EXP(UULMXR)
                  IF(UUMIN.LE.0.) THEN
                     UULMNR=-100
                  ELSE
                     UULMNR = (LOG(UUMIN)-TVLEFF)*(NEXTOT-1)
                  ENDIF
                  IF(UULMNR.LT.-88) THEN
                     UUMNR = 0.
                  ELSE
                     UUMNR = EXP(UULMNR)
                  ENDIF
               ENDIF
               AITOT  = NPART (JEMISS) * FLAMMX * ( UUMXR - UUMNR )
               AIPRO  = NPART (JEMISS) * FLAMMX * MIN ( UUMXR, ONEONE )
 2800       CONTINUE
            NTENT  = INT  (AIPRO)
            CALL GRNDM(RNDM,1)
            RNTENT = RNDM (1)
            IF ( RNTENT .LT. AIPRO - NTENT ) NTENT = NTENT + 1
 3000       CONTINUE
            ITACC = 0
            DO 4100 IT = 1, NTENT
               GO TO ( 3100, 3200, 3300, 3400, 3500, 3600, 3700 ),
     &            ISAMPL
               GO TO 4000
 3100          CONTINUE
                  CALL GRNDM(RNDM,1)
                  UURND = AEFRPA * ( 1.D+00 - AITOT / FLAMMX
     &                  * RNDM (1) )**0.3333333333333333D+00
                  UURND = AEFRPA - UURND + UUMIN
               GO TO 3800
 3200          CONTINUE
                  CALL GRNDM(RNDM,1)
                  RNDEPS = RNDM (1)
                  IF ( RNDEPS .LT. AITO1 ) THEN
                     RNDEPS = RNDEPS / AITO1
                     UURND  = SQRT ( RNDEPS * ( UUDV2 - UUMN2 ) + UUMN2)
                  ELSE
                     RNDEPS = RNDEPS - AITO1
                     RNDEPS = AIHLP + RNDEPS * AITOT * EDENOM / FLAMMX
                     UURND  = 2.D+00 * UUDIV - SQRT ( UUDIV**2
     &                      - 2.D+00 * RNDEPS )
                  END IF
               GO TO 3800
 3300          CONTINUE
                  CALL GRNDM(RNDM,1)
                  RNDEPS = RNDM (1)
                  IF ( RNDEPS .LT. AITO1 ) THEN
                     RNDEPS = RNDEPS * AITOT
                     UURND  = RNDEPS / FCHLP + UUMIN**2
                     UURND  = SQRT (UURND)
                  ELSE
                     RNDEPS = ( RNDEPS - AITO1 ) * AITOT + AIHLP
                     UURND  = RNDEPS / FCHLP / UUDIV + UUDIV
                  END IF
               GO TO 3800
 3400          CONTINUE
                  CALL GRNDM(RNDM,1)
                  UURND = TVEFF * ( RNDM (1) * AITOT / FLAMMX
     &                  / FEMAX / NPART(JEMISS) + UUMNR )**
     &                  ( 1.D+00 / (NEXTOT-1) )
                  IF ( UURND .LE. AEFRPA ) THEN
                     LSPREJ = .FALSE.
                     FREJE  = 1.D+00
                  ELSE IF ( UURND .LE. 2.D+00 * AEFRPA ) THEN
                     FREJE  = ( UURND**2 - 2.D+00 * ( UURND - AEFRPA )
     &                      **2 ) / UURND**2
                  ELSE
                     FREJE  = 2.D+00 * AEFRPA**2 / UURND**2
                  END IF
               GO TO 3800
 3500          CONTINUE
                  CALL GRNDM(RNDM,1)
                  RNDEPS = RNDM (1)
                  IF ( RNDEPS .LT. AITO1 ) THEN
                     RNDEPS = RNDEPS / AITO1
                     UURND  = RNDEPS * ( UUDV3 - UUMN3 ) + UUMN3
                     UURND  = UURND**0.333333333333333D+00
                     LSPREJ = .FALSE.
                  ELSE IF ( RNDEPS .LT. AITO1 + AITO2 ) THEN
                     RNDEPS = RNDEPS - AITO1
                     RNDEPS = RNDEPS * AITOT * EDENOM / FLAMMX + AIHLP
                     UURND  = 0.5D+00 * ( UUDIV + SQRT ( UUDIV**2
     &                      + 1.333333333333333D+00 * RNDEPS / FCHLP
     &                      / UUDIV ) )
                     LSPREJ = .FALSE.
                  END IF
               GO TO 3800
 3600          CONTINUE
                  CALL GRNDM(RNDM,1)
                  RNDEPS = RNDM (1)
                  IF ( RNDEPS .LT. AITO1 ) THEN
                     RNDEPS = RNDEPS / AITO1
                     UURND  = RNDEPS * ( UUDV4 - UUMN4 ) + UUMN4
                     UURND  = UURND**0.25D+00
                     LSPREJ = .FALSE.
                  ELSE IF ( RNDEPS .LT. AITO1 + AITO2 ) THEN
                     RNDEPS = RNDEPS - AITO1
                     RNDEPS = RNDEPS * AITOT * EDENOM / FLAMMX + AIHLP
                     UURND  = RNDEPS / FCHLP + UUDV4
                     UURND  = UURND**0.25D+00
                     LSPREJ = .TRUE.
                     FREJE  = 1.D+00 - 2.D+00 * ( 1.D+00 - UUDIV / UURND
     &                      )**3
                  ELSE
                     RNDEPS = RNDEPS - AITO1 - AITO2
                     RNDEPS = RNDEPS * AITOT * EDENOM / FLAMMX + AIHL2
                     UURND  = AEFRPA + SQRT ( AEFRPA**2 + RNDEPS /
     &                      ( 12.D+00 * FCHLP * AEFRPA**2 ) )
                     LSPREJ = .FALSE.
                  END IF
               GO TO 3800
 3700          CONTINUE
                  CALL GRNDM(RNDM,1)
                  UURND = TVEFF * ( RNDM (1) * AITOT / FLAMMX
     &                  / NPART(JEMISS) + UUMNR )**
     &                  ( 1.D+00 / (NEXTOT-1) )
 3800          CONTINUE
               IF ( LSPREJ ) THEN
                  CALL GRNDM(RNDM,1)
                  IF ( RNDM (1) .GE. FREJE ) GO TO 4100
               END IF
               EPS    = TVEFF - UURND - BNDGAV
 4000          CONTINUE
               IF ( LPRFIX ) THEN
                  ITACC = ITACC + 1
                  EPSFIX (ITACC) = EPS - EPSMIN + BNDGEN
               END IF
               DEPSEX = EPS - EPSEEX (JEMISS)
               IF ( DEPSEX .GT. 0.D+00 ) THEN
                  IF ( EEXDLT (JEMISS) .LT. EEXMNM (JEMISS) .AND.
     &               DEPSEX .GT. 0.5D+00 * ( EPSDLT (JEMISS) - EPSEEX
     &              (JEMISS) ) ) THEN
                     DEPSEX = EPS - EPSDLT (JEMISS)
                     IF ( DEPSEX .GT. 0.5D+00 * ( EPSDLT (JEMISS)
     &                  - EPSANY (JEMISS) ) ) THEN
                        EPS = EPSANY (JEMISS)
                     ELSE
                        EPS = EPSDLT (JEMISS)
                     END IF
                  ELSE
                     EPS = EPSEEX (JEMISS)
                  END IF
                  IF ( EPS .LE. EPSMIN ) GO TO 4100
               END IF
               AMNHLP = UMO - EPS + EPSMIN - AMNUCL (IEMISS)
               ERNCM (JEMISS) = 0.5D+00 * ( UMO2 + AMNHLP**2
     &                        - AMNUSQ (IEMISS) ) / UMO
               EPS    = UMO - ERNCM (JEMISS) - AMNUCL (IEMISS)
               EPSINV = 0.5D+00 * ( UMO2 - ( AMNUCL (IEMISS)
     &                + AMNHLP )**2 ) / AMNHLP
               EPSWLL = EPSINV + AEFRAV + BNDUSE + AMNUCL (IEMISS)
               BETWLL = SQRT ( 1.D+00 - AMNUSQ (IEMISS) / EPSWLL**2 )
               EKEWLL = EPSWLL - AMNUCL (IEMISS)
               EPSCOU = UMO - AMNHLP - AMNUCL (IEMISS)
               IF ( JEMISS .EQ. 1 ) THEN
                  GSNGLV = 1.5D+00 * ( ANOW - ZNOW ) * SQRT ( ( EKEWLL
     &                   + EKFCEN (IEMISS) - AEFRAV ) / EKFCEN (IEMISS))
     &                   / EKFCEN (IEMISS)
                  AALPHA = 0.76D+00 + 2.2D+00 / RMASS (IBRES-1)
                  BBETAA = ( 2.22D+00 / RMASS (IBRES-1)**2 - 0.05D+00 )
     &                / AALPHA * 1.D-03
                  SIGINV = 0.1D+00 * XINNEU ( EPSINV, ZAFTR (JEMISS),
     &                                        BBETAA )
                  IF ( EKEWLL .LE. 0.04D+00 ) THEN
                     SIGMAP = 3.D+03 * PI / ( 1.206D+03 * EKEWLL + (
     &                      -1.86D+00 + 0.09415D+03 * EKEWLL
     &                      + 0.0001306D+06 * EKEWLL**2 )**2 ) + 1.D+03
     &                      * PI / ( 1.206D+03 * EKEWLL + ( 0.4223D+00
     &                      + 0.13D+03 * EKEWLL )**2 )
                     IF ( EKEWLL .LT. 0.02D+00 ) THEN
                        SIGMAN = 0.3333333333333333D+00 * SIGMAP
                     ELSE
                        SIGMAN = 10.63D+00 / BETWLL**2 - 29.92D+00
     &                         / BETWLL + 42.9D+00
                     END IF
                  ELSE
                     SIGMAP = 34.10D+00 / BETWLL**2 - 82.2D+00 / BETWLL
     &                      + 82.2D+00
                     SIGMAN = 10.63D+00 / BETWLL**2 - 29.92D+00 / BETWLL
     &                      + 42.9D+00
                  END IF
                  SIGMAP = 0.1D+00 * SIGMAP
                  SIGMAN = 0.1D+00 * SIGMAN
                  SGNUNU = ( ( ANOW - ZNOW - ANNEUT ) * SIGMAN
     &                   + ( ZNOW - ANPROT ) * SIGMAP )
     &                   / ( ANOW - ANPTOT )
               ELSE
                  GSNGLV = 1.5D+00 * ZNOW * SQRT ( ( EKEWLL
     &                   + EKFCEN (IEMISS) - AEFRAV ) / EKFCEN (IEMISS))
     &                   / EKFCEN (IEMISS)
                  SIGINV = 1.D-01 * XINPRO ( EPSINV, ZAFTR (JEMISS),
     &                                       ETHRES )
                  IF ( EKEWLL .LE. 0.04D+00 ) THEN
                     SIGMAN = 3.D+03 * PI / ( 1.206D+03 * EKEWLL + (
     &                      -1.86D+00 + 0.09415D+03 * EKEWLL
     &                      + 0.0001306D+06 * EKEWLL**2 )**2 ) + 1.D+03
     &                      * PI / ( 1.206D+03 * EKEWLL + ( 0.4223D+00
     &                      + 0.13D+03 * EKEWLL )**2 )
                     IF ( EKEWLL .LT. 0.02D+00 ) THEN
                        SIGMAP = 0.3333333333333333D+00 * SIGMAN
                     ELSE
                        SIGMAP = 10.63D+00 / BETWLL**2 - 29.92D+00
     &                         / BETWLL + 42.9D+00
                     END IF
                  ELSE
                     SIGMAN = 34.10D+00 / BETWLL**2 - 82.2D+00 / BETWLL
     &                      + 82.2D+00
                     SIGMAP = 10.63D+00 / BETWLL**2 - 29.92D+00 / BETWLL
     &                      + 42.9D+00
                  END IF
                  SIGMAP = 0.1D+00 * SIGMAP
                  SIGMAN = 0.1D+00 * SIGMAN
                  SGNUNU = ( ( ANOW - ZNOW - ANNEUT ) * SIGMAN
     &                   + ( ZNOW - ANPROT ) * SIGMAP )
     &                   / ( ANOW - ANPTOT )
               END IF
               ZITA   = AEFRAV / EKEWLL
               IF ( ZITA .LE. 0.5D+00 ) THEN
                  PZITA = 1.D+00 - 1.4D+00 * ZITA
               ELSE
                  PZITA = 1.D+00 - 1.4D+00 * ZITA + 0.4D+00 * ZITA
     &                  * ( 2.D+00 - 1.D+00 / ZITA )**2.5D+00
               END IF
               RATEC  = SIGINV * RATEC0 * EPSINV / GSNGLV
               IF ( RATEC .LE. 0.D+00 ) GO TO 4100
               ALAMDC = BETWLL / RATEC
               ALAMDI = 1.D+00 / ( SGNUNU * RHONUC * PZITA ) + MAX (
     &                  TWOTWO * PI * PLABRC * BETWLL * EPSWLL
     &                / AMNUSQ (IEMISS), RZNUCL )
               AMUTOT = 1.D+00 / ALAMDI + 1.D+00 / ALAMDC
               RATEI  = BETWLL / ALAMDI
               FLAMDA = RATEC / ( RATEC + CLAMDI * RATEI ) / FLAMMX
               CALL GRNDM(RNDM,1)
               IF ( RNDM (1) .LT. FLAMDA ) GO TO 4200
 4100       CONTINUE
            IF ( JEMISS .EQ. JEMIS1 ) GO TO 4500
            IF ( ITACC .LE. NPART (JEMISS) ) THEN
               NPRFIX (JEMISS) = NPRFIX (JEMISS) + ITACC
               DO 4150 IT = 1, ITACC
                  EPRFIX = EPRFIX + EPSFIX (IT)
 4150          CONTINUE
            ELSE
               NPRFIX (JEMISS) = NPART (JEMISS)
               ACCEP = 0.D+00
               DO 4160 IT = 1, ITACC
                  PRACC  = ( NPART (JEMISS) - ACCEP )
     &                   / ( ITACC - IT + 1 )
                  CALL GRNDM(RNDM,1)
                  IF ( RNDM (1) .LT. PRACC ) THEN
                     EPRFIX = EPRFIX + EPSFIX (IT)
                     ACCEP  = ACCEP  + 1.D+00
                  END IF
 4160          CONTINUE
            END IF
            GO TO 4500
 4200       CONTINUE
            ZNOW = ZAFTR (JEMISS)
            ANOW = AAFTR
            IBRES = IBRES - 1
            ICRES = ICRES - JEMISS + 1
            AMMRES = AMMAFT (JEMISS)
            AMNRES = AMNAFT (JEMISS)
            EEPCM  = EPS + AMNUCL (IEMISS)
            ERESCM = UMO - EEPCM
            IF ( LNUCTS ) THEN
            ELSE IF ( ICYCL .EQ. 1 .AND. LEMISS ) THEN
               C (1) = C2ND (1)
               C (2) = C2ND (2)
               C (3) = C2ND (3)
            ELSE IF ( ICYCL .LE. 1 .AND. PCMORI .LE. ANGLGB ) THEN
               C (1) = C1ST (1)
               C (2) = C1ST (2)
               C (3) = C1ST (3)
            ELSE IF ( ICYCL .LE. 1 ) THEN
               PCMSXT = - AMNUCL (IEMISS) * ETAX
               PCMSYT = - AMNUCL (IEMISS) * ETAY
               PCMSZT = - AMNUCL (IEMISS) * ETAZ
               APSHAR = MAX ( ANPTOT + NEMISS - 1, ONEONE )
               PCMSX  = ( CXAXCM * PCMORI + PCMSXT ) / APSHAR - PCMSXT
               PCMSY  = ( CYAXCM * PCMORI + PCMSYT ) / APSHAR - PCMSYT
               PCMSZ  = ( CZAXCM * PCMORI + PCMSZT ) / APSHAR - PCMSZT
*
               PCMSTT = PCMSXT**2 + PCMSYT**2 + PCMSZT**2
               PCMSPT = PCMSX**2  + PCMSY**2  + PCMSZ**2
               ECMSTT = SQRT ( AMNUSQ (IEMISS) + PCMSTT )
               ECMSPT = SQRT ( AMNUSQ (IEMISS) + PCMSPT )
*
               UMMO2  = 2.D+00 * AMNUSQ (IEMISS) + 2.D+00 * ECMSTT
     &                * ECMSPT - 2.D+00 * ( PCMSXT * PCMSX + PCMSYT
     &                * PCMSY + PCMSZT * PCMSZ )
               UMMO   = SQRT ( UMMO2 )
               GAMCCM = ( ECMSTT + ECMSPT ) / UMMO
               ETAXM  = ( PCMSX + PCMSXT ) / UMMO
               ETAYM  = ( PCMSY + PCMSYT ) / UMMO
               ETAZM  = ( PCMSZ + PCMSZT ) / UMMO
               ETAPCM = ETAXM * PCMSX + ETAYM * PCMSY + ETAZM * PCMSZ
               PHELP  = ECMSPT - ETAPCM / ( GAMCCM + 1.D+00 )
               PCCMSX = PCMSX - ETAXM * PHELP
               PCCMSY = PCMSY - ETAYM * PHELP
               PCCMSZ = PCMSZ - ETAZM * PHELP
               PCCMSP = SQRT ( PCCMSX**2 + PCCMSY**2 + PCCMSZ**2 )
               CALL RACO (C(1),C(2),C(3))
               SCAPRO = PCCMSX * C (1) + PCCMSY * C (2) + PCCMSZ * C (3)
               IF ( SCAPRO .LT. 0.D+00 ) THEN
                  C (1) = - C (1)
                  C (2) = - C (2)
                  C (3) = - C (3)
               END IF
               PCCMSX = C(1) * PCCMSP
               PCCMSY = C(2) * PCCMSP
               PCCMSZ = C(3) * PCCMSP
               ETAPCM = ETAXM * PCCMSX + ETAYM * PCCMSY + ETAZM * PCCMSZ
               PHELP  = 0.5D+00 * UMMO + ETAPCM / ( GAMCCM + 1.D+00 )
               PCMSX  = PCCMSX + ETAXM * PHELP
               PCMSY  = PCCMSY + ETAYM * PHELP
               PCMSZ  = PCCMSZ + ETAZM * PHELP
               PCMSPT = SQRT ( PCMSX**2 + PCMSY**2 + PCMSZ**2 )
               C (1)  = PCMSX / PCMSPT
               C (2)  = PCMSY / PCMSPT
               C (3)  = PCMSZ / PCMSPT
               ETAPCM = - ETAPCM
               PHELP  = 0.5D+00 * UMMO + ETAPCM / ( GAMCCM + 1.D+00 )
               PCMSX  = - PCCMSX + ETAXM * PHELP
               PCMSY  = - PCCMSY + ETAYM * PHELP
               PCMSZ  = - PCCMSZ + ETAZM * PHELP
               PCMSPT = SQRT ( PCMSX**2 + PCMSY**2 + PCMSZ**2 )
               C2ND (1) = PCMSX / PCMSPT
               C2ND (2) = PCMSY / PCMSPT
               C2ND (3) = PCMSZ / PCMSPT
            ELSE
               UMULEG (0) = 1.0D+00
               UMULEG (1) = 0.6666666666666667D+00
               UMULEG (2) = 0.25D+00
               ACOLEG (0) = 0.5D+00
               ACOLEG (1) = 1.5D+00 * UMULEG (1)**ICYCL
               ACOLEG (2) = 2.5D+00 * UMULEG (2)**ICYCL
               COSTHE = COSLEG ( ACOLEG )
               SINTHE = SQRT ( ( 1.D+00 - COSTHE ) * ( 1.D+00 + COSTHE )
     &                )
 4300          CONTINUE
                  CALL GRNDM(RNDM,2)
                  RPHI1 = 2.D+00 * RNDM (1) - 1.D+00
                  RPHI2 = 2.D+00 * RNDM (2) - 1.D+00
                  RPHI12 = RPHI1 * RPHI1
                  RPHI22 = RPHI2 * RPHI2
                  RSQ = RPHI12 + RPHI22
               IF ( RSQ .GT. 1.D+00 ) GO TO 4300
               SINPHI = 2.D+00 * RPHI1 * RPHI2 / RSQ
               COSPHI = ( RPHI12 - RPHI22 ) / RSQ
               SINT02 = ( 1.D+00 - CZAXCM ) * ( 1.D+00 + CZAXCM )
               IF ( SINT02 .LT. ANGLSQ ) THEN
                  C (1) = COSPHI * SINTHE
                  C (2) = SINPHI * SINTHE
                  C (3) = CZAXCM * COSTHE
               ELSE
                  SINTH0 = SQRT ( SINT02 )
                  UPRIME = SINTHE * COSPHI
                  VPRIME = SINTHE * SINPHI
                  COSPH0 = CXAXCM / SINTH0
                  SINPH0 = CYAXCM / SINTH0
                  C (1) = UPRIME * COSPH0 * CZAXCM - VPRIME * SINPH0
     &                  + COSTHE * CXAXCM
                  C (2) = UPRIME * SINPH0 * CZAXCM + VPRIME * COSPH0
     &                  + COSTHE * CYAXCM
                  C (3) = - UPRIME * SINTH0 + COSTHE * CZAXCM
               END IF
            END IF
            PCMS  = SQRT ( ( EEPCM - AMNUCL (IEMISS) ) * ( EEPCM
     &            + AMNUCL (IEMISS) ) )
            UMO2   = ( ERESCM - PCMS ) * ( ERESCM + PCMS )
            UMO    = SQRT (UMO2)
            PCMSX = PCMS * C (1)
            PCMSY = PCMS * C (2)
            PCMSZ = PCMS * C (3)
            ETAPCM = PCMSX * ETAX + PCMSY * ETAY + PCMSZ * ETAZ
            NP = NP + 1
            KPART (NP) = 1 - 7 * ( JEMISS - 2 )
            TKI   (NP) = GAMCM * EEPCM + ETAPCM - AMNUCL (IEMISS)
            PHELP = ETAPCM / ( GAMCM + 1.D+00 ) + EEPCM
            WEI   (NP) = WEE
            PLBPX = PCMSX + ETAX * PHELP
            PLBPY = PCMSY + ETAY * PHELP
            PLBPZ = PCMSZ + ETAZ * PHELP
            PHELP = SQRT ( PLBPX * PLBPX + PLBPY * PLBPY
     &            + PLBPZ * PLBPZ )
            CXR   (NP) = PLBPX / PHELP
            CYR   (NP) = PLBPY / PHELP
            CZR   (NP) = PLBPZ / PHELP
            PLR   (NP) = PHELP
            IGREYP = IGREYP + JEMISS - 1
            IGREYN = IGREYN + 2 - JEMISS
            IF ( JEMISS .EQ. 1 ) THEN
               ISTORE = IGREYN
            ELSE
               ISTORE = IGREYP
            END IF
            IF ( LNUCTS ) THEN
               NPCYCL (ISTORE,JEMISS) = -ICYCL
            ELSE
               NPCYCL (ISTORE,JEMISS) = ICYCL
            END IF
            EINTR  = EINTR  + TKI (NP) + AMNUCL (IEMISS)
            PXINTR = PXINTR + CXR (NP) * PLR (NP)
            PYINTR = PYINTR + CYR (NP) * PLR (NP)
            PZINTR = PZINTR + CZR (NP) * PLR (NP)
            IBINTR = IBINTR + IBAR ( KPART (NP) )
            ICINTR = ICINTR + ICH  ( KPART (NP) )
            ERES  = GAMCM * ERESCM - ETAPCM
            EKRES = ERES - UMO
            PHELP = - ETAPCM / ( GAMCM + 1.D+00 ) + ERESCM
            PXRES = - PCMSX + ETAX * PHELP
            PYRES = - PCMSY + ETAY * PHELP
            PZRES = - PCMSZ + ETAZ * PHELP
            PTRES2= PXRES * PXRES + PYRES * PYRES + PZRES * PZRES
            PTRES = SQRT (PTRES2)
            TVCMS = UMO - AMNRES
            IF ( TVCMS .LT. 0.D+00 ) THEN
               IF ( TVCMS .LT. - GAMMIN ) THEN
                  WRITE (LUNOUT,*)' PREPRE: TVCMS',TVCMS
                  WRITE (LUNERR,*)' PREPRE: TVCMS',TVCMS
                  WRITE (77,*)' ^^^ PREPRE: TVCMS',TVCMS
               END IF
               TVCMS = 0.D+00
            END IF
            EHLFIX = MIN ( TVCMS, EHLFIX )
            TVEFF = TVCMS - EHLFIX
            IF ( TVEFF .LE. 0.D+00 ) GO TO 7000
            GAMCM = ERES  / UMO
            ETACM = PTRES / UMO
            ETAX  = PXRES / UMO
            ETAY  = PYRES / UMO
            ETAZ  = PZRES / UMO
            ETAPCM = ETAX * PXORI + ETAY * PYORI + ETAZ * PZORI
            PHELP = EKORI + AM (KPORI) - ETAPCM / ( GAMCM + 1.D+00 )
            PCMSX = PXORI - ETAX * PHELP
            PCMSY = PYORI - ETAY * PHELP
            PCMSZ = PZORI - ETAZ * PHELP
            PCMORI = PCMSX**2 + PCMSY**2 + PCMSZ**2
            IF ( PCMORI .LE. ANGLGB ) THEN
               PCMORI = 0.D+00
               CALL RACO ( CXAXCM, CYAXCM, CZAXCM )
            ELSE
               PCMORI = SQRT ( PCMORI )
               CXAXCM = PCMSX / PCMORI
               CYAXCM = PCMSY / PCMORI
               CZAXCM = PCMSZ / PCMORI
            END IF
            IAAFT = IBRES-1
            IZAFT = ICRES
            AAFTR = IAAFT
            ZAFTR (1) = IZAFT
            AMMAFT(1) = AAFTR * AMUAMU + 1.D-03
     &                * FKENER (AAFTR,ZAFTR(1))
            AMNAFT(1) = AMMAFT (1) - ZAFTR (1) * AMELEC
     &                + ELBNDE (IZAFT)
            CALL EEXLVL ( IAAFT, IZAFT, EEXDLT (1), EEXMNM (1), EEXDUM )
            BNDGAV = BNENRG (2)
            EPSMX  (1) = UMO - AMNRES - BNDGAV - EHLFIX
            EPSDLT (1) = UMO - AMNRES - BNDGAV - EEXDLT (1)
            EPSEEX (1) = UMO - AMNRES - BNDGAV - EEXMNM (1)
            IZAFT = ICRES - 1
            ZAFTR (2) = IZAFT
            AMMAFT(2) = AAFTR * AMUAMU + 1.D-03
     &                * FKENER (AAFTR,ZAFTR(2))
            AMNAFT(2) = AMMAFT (2) - ZAFTR (2) * AMELEC
     &                + ELBNDE (IZAFT)
            CALL EEXLVL ( IAAFT, IZAFT, EEXDLT (2), EEXMNM (2), EEXDUM )
            BNDGAV = BNENRG (1)
            EPSMX  (2) = UMO - AMNRES - BNDGAV - EHLFIX
            EPSDLT (2) = UMO - AMNRES - BNDGAV - EEXDLT (2)
            EPSEEX (2) = UMO - AMNRES - BNDGAV - EEXMNM (2)
            EEXANW (1) = 0.D+00
            EPSANY (1) = EPSMX (1)
            EEXANW (2) = 0.D+00
            EPSANY (2) = EPSMX (2)
            NPART (JEMISS) = NPART (JEMISS) - 1
            IF ( JEMFIX .GT. 0 ) THEN
               NPART (1) = NPART (1) + NPRFIX (1)
               NPART (2) = NPART (2) + NPRFIX (2)
            END IF
            NEMISS = NEMISS + 1
            LEMISS = .TRUE.
            TEMNU = SQRT ( TVCMS / ANOW * ALEVEL  )
            IF ( EPSMX (1) + EPSMX (2) .LT. 2.D+00 * TEMNU )
     &         GO TO 7000
            SIGIN0 = PI * ( R0SIGK * RMASS (IBRES-1) )**2
            NEXMX = SQRT ( 2.D+00 * ANOW * TVCMS / ALEVEL )
            ANPROT =  NPART (2)
            ANNEUT =  NPART (1)
            NPTOT  = NPART (1) + NPART (2)
            ANPTOT = NPTOT
            GO TO 6000
 4500    CONTINUE
 4600    CONTINUE
         IF ( LNUCTS ) THEN
            LEMISS = .TRUE.
            GO TO 6000
         ELSE
            LEMISS = .FALSE.
            NPART (1) = NPART (1) + NPRFIX (1)
            NPART (2) = NPART (2) + NPRFIX (2)
            IF ( JEMFIX .GT. 0 ) THEN
               TVEFF = TVCMS - EHLFIX
               EPSMX (JEMFIX) = EPSMX0
               ERNCM (JEMFIX) = ERNCM0
            END IF
         END IF
 5000    CONTINUE
         ANPROT =  NPART (2)
         ANNEUT =  NPART (1)
         IF ( NPTOT .GT. 0 ) THEN
            PNPROT = ( ZNOW - ANPROT ) * ( 3.D+00 * ANNEUT + ANPROT )
     &             / ( ANPROT * ( ZNOW - ANPROT + 3.D+00 * ( ANOW
     &             - ANNEUT - ZNOW ) ) + ANNEUT * ( 3.D+00 * ( ZNOW
     &             - ANPROT ) + ANOW - ANNEUT - ZNOW ) )
         ELSE
            PNPROT = ZNOW / ANOW
         END IF
         CALL GRNDM(RNDM,1)
         IF ( RNDM (1) .LT. PNPROT ) THEN
            NPART (2) = NPART (2) + 1
            ZCOLL = ZCOLL - 1.D+00
         ELSE
            NPART (1) = NPART (1) + 1
         END IF
         NPTOT  = NPART (1) + NPART (2)
         ANPTOT = NPTOT
         NHOLE  = NHOLE + 1
         ACOLL  = ACOLL - 1.D+00
         AVEBIN = ( ( BBTAR - ZZTAR - ACOLL + ZCOLL ) * BNENRG (2)
     &          + ( ZZTAR - ZCOLL ) * BNENRG (1) ) / ( BBTAR - ACOLL )
 6000    CONTINUE
         IF ( LNUCTS ) THEN
            RHOIMP = RHOSAV
            EKFIMP = EKFSAV
         END IF
      IF ( NHOLE + NHLFIX + NPTOT .LT. NEXMX .AND.
     &     NPTOT .LE. NINT (0.5D+00 * ANOW) ) GO TO 1000
 7000 CONTINUE
*=== End of subroutine prepre =========================================*
      RETURN
      END
#if defined(CERNLIB_HPUX)
$OPTIMIZE ON
#endif
